function keffmg1_5_mod_with_hebert
disp('----------------------------------------')
disp('------- 1D mesh edge diffusion ---------')
disp('----------------------------------------')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  multi-group diffusion solver for slab geometry,  j.roberts            %
%  this solver uses mesh edge differencing                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BOUNDARY CODITIONS                                                     %
%   0 == vacuum condition (boundary flux = 0)                            %
%   1 == reflective condition                                            %
%   2 == zero incoming current using partial current (better than 0)     %
% SOURCE OPTIONS                                                         %
%   0 == fixed-source (volume source)                                    %
%   1 == eigenvalue  (solves via power iteration)                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; %clc;

bcL = 2;        % left boundary condition (do not mix 0 with 1 or 2!)
bcR = 2;        % right boundary condition
it = 1;         % source option

% material and cross-section parameters
% cross-section format:
%    for all materials
%        for all groups
%           read difc(m,g) siga(m,g) nusi(m,g) sigsg1->gi sigsg2->gi ...
% we limit the data only to downscattering

% Simple 1-D problem
numg = 2; % number of groups
numm = 2; % number of materials
file = 'xsec_2g';
[difc,siga,nusi,sigs] = rdxsec(file,numm,numg);

% mesh information
coarse  =   [ 0   20];
fine    =   [   1   ]*50;
% material placement
matl    =   [   1   ];
% source in each coarse mesh
source    = [   0         
                0   ]';       
          
numc   = length(coarse)-1;   % number of course meshes
numf   = sum(fine);          % total number of fine meshes

% initialize the finemesh data values
dd = zeros(numf,numg);      % diffusion coefficient
aa = zeros(numf,numg);      % absorption cross-section
nf = zeros(numf,numg);      % nu*fission cross-section
sc = zeros(numf,numg,numg); % 
ss = zeros(numf,numg);
dx = zeros(numf,1);

j = 0;
for i = 1:numc
    dx( (j+1):(j+fine(i)) ) = (coarse(i+1)-coarse(i))/fine(i);
    for g = 1:numg
	    dd( (j+1):(j+fine(i)), g ) = difc( matl(i), g  ); 
	    aa( (j+1):(j+fine(i)), g ) = siga( matl(i), g  );    
	    nf( (j+1):(j+fine(i)), g ) = nusi( matl(i), g  ); 
        for k = (j+1):(j+fine(i))
            sc( k, g, : ) = sigs( matl(i), g, :  );
        end
        if it == 0
            ss( (j+1):(j+fine(i)), g  ) = source(i, g ); 
        end
    end
    j = sum(fine(1:i));
end

if sum( dx >= sqrt(dd(:,end)./aa(:,end)) ) > 0
    disp('meshes are too large!') % want dx < L = sqrt(DiffCoef/SigA)
end

if bcL == 0 % then bcR must be 0, too. 
    n = numf-1; % (# MESHES)-1 flux points
elseif bcL == 1 || bcL == 2
    n = numf+1; % (# MESHES)+1 flux points
else
    disp('*** bc unrecognized or is not yet supported ***')
end
lala=1;
%%%%%%%%%%%%%%%%%%%% COEFFICIENT MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% build the coefficient matrix (just diagonals)
AD = zeros(n,numg);
AL = zeros(n-1,numg);
AU = zeros(n-1,numg);
for g = 1:numg
    if bcL == 0  % zero boundary flux (no extrapolation)
        ii=0;
        AU(1,g)   = -dd(2,g)/dx(2,1);    
        AD(1,g)   = -AU(1,g) + dd(1,g)/dx(1,1) + ...
                     0.25*(aa(2,g)*dx(2,1)+aa(1,g)*dx(1,1));
    elseif bcL == 1 % reflective
        ii=1;
        AU(1,g)   = -dd(1,g)/dx(1,1);
        AD(1,g)   = -AU(1,g) + 0.5*(aa(1,g)*dx(1,1)) ;
    else % better vacuum approximation
        ii=1;
        AU(1,g)   = -dd(1,g)/dx(1,1);
        AD(1,g)   = -AU(1,g) + 0.5*(1+lala*dx(1)*aa(1,g));
    end
    if bcR == 0  % zero boundary flux (no extrapolation)
        ii=0;
        AL(n-1,g) = -dd(n,g)/dx(n,1); 
        AD(n,g)   = -AL(n-1,g) + dd(n+1,g)/dx(n+1,1) + ...
                     0.25*(aa(n+1,g)*dx(n+1,1)+aa(n,g)*dx(n,1));
    elseif bcR == 1 % reflective
        ii=1;
        AL(n-1,g) = -dd(n-1,g)/dx(n-1,1);
        AD(n,g)   = -AL(n-1,g) + ...
                    0.5*(aa(n-1,g)*dx(n-1,1));
    else % better vacuum approximation
        ii=1;
        AL(n-1,g) = -dd(n-1,g)/dx(n-1,1);
        AD(n,g)   = -AL(n-1,g) + 0.5*(1+lala*aa(n-1,g)*dx(n-1,1));
    end  
    for i = (2-ii):(n-1-ii)  % flux "1" to "100"
        AL(i+ii-1,g) = -dd(i,g)/dx(i,1);
        AU(i+ii,g)   = -dd(i+1,g)/dx(i+1,1);    
        AD(i+ii,g)   = -( AL(i+ii-1,g) + AU(i+ii,g) ) + ...
            0.5*(aa(i+1,g)*dx(i+1,1)+aa(i,g)*dx(i,1));
    end
end

if it == 1
    nusig = zeros(n,numg);
else
    s = zeros(n,1);
end
dxx = zeros(n,1);
sct = zeros(n,numg,numg);


%%%%%%%%%%%%%%%%%%  CONSTANTS FOR BOUNDARIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for g = 1:numg
    if bcL > 0 % set values for constants at boundaries
        if it == 1
            if bcL == 1
                nusig(1,g) = 0.5*dx(1)*nf(1,g);
            else
                nusig(1,g) = lala*0.5*dx(1)*nf(1,g);
            end
            if bcR == 1
                nusig(n,g) = 0.5*dx(n-1)*nf(n-1,g);  
            else
                nusig(n,g) = lala*0.5*dx(n-1)*nf(n-1,g);
            end
        else    
            if bcL == 1
                s(1,g) = 0.5*dx(1)*ss(1,g);
            else
                s(1,g) = 0.0;
            end
            if bcR == 1
                s(n,g) = 0.5*dx(n-1)*ss(n-1,g);
            else
                s(n,g) = 0.0;
            end
        end
        dxx(1)   = 0.5*dx(1);
        dxx(n)   = 0.5*dx(n-1);
        sct(1,g,:)   = 0.5*dx(1)*sc(1,g,:);
        sct(n,g,:)   = 0.5*dx(n-1)*sc(n-1,g,:);
    end
    for i = (1):(n-ii*2) 
        if it == 1
            nusig(i+ii,g) = 0.5*dx(i)*nf(i+1,g)+0.5*dx(i+1)*nf(i+1,g);
        else
            s(i+ii,g)     = 0.5*dx(i)*ss(i,g)+0.5*dx(i+1)*ss(i+1,g);
        end
        sct(i+ii,g,:)   = 0.5*dx(i)*sc(i+1,g,:) + 0.5*dx(i+1)*sc(i,g,:);
    end
end
for i = 1:(n-ii*2)
    dxx(i+ii,:) = 0.5*(dx(i+1,:)+dx(i,:));
end

if it == 1
    sindx = find( sum(nusig(:,:)') ); % if no nusig, so source  
    s = 0.4*ones(n,1).*( sum(nusig(:,:)' > 0 ) )';  
    s = 0.01 * s / max(s);
end

%------------- make the x-vector for plotting------------------------
x = zeros(length(dx)+1,1);
for i = 2:1:length(dx)+1
	x(i) = x(i-1,1) + dx(i-1,1);
end

phi     = zeros(n,g);
scsrc   = zeros(n,1);

if it == 1  % eigenvalue search -------------------------------------
    [k s phi] = keffit( AU, AL, AD, dxx, nusig, sct, s, 1e-8, 1e-9,100);
else        % fixed source, probably want function ------------------
    for g = 1:numg
        % compute scattering source
        for k = 1:n
            for gg = 1:g-1 % only down scatter
                scsrc(k,1) = scsrc(k,1) + sct(k,g,gg)*phi(k,gg);
            end
        end
        % solve for phi_g
        phi(:,g) = tridiag(AU(:,g),AL(:,g),AD(:,g),s(:,g)+scsrc(:,1));
    end  
end


figure(1)
hold on
if numg==2
    if bcL == 0 && bcR == 0
        plot(x,[0 phi(:,1)' 0],'k',x,[0 phi(:,2)' 0],'b','LineWidth',2)
        axis([min(x) max(x) 0 1.25*max(max(phi))]);
        xlabel('{x(cm)}')
        ylabel('\phi(x)')
    else
        plot(x,phi(:,1),'k',x,phi(:,2),'b','LineWidth',2)
        axis([min(x) max(x) 0 1.25*max(max(phi))]);
        xlabel('{x(cm)}')
        ylabel('\phi(x)')
    end
elseif numg==1
    if bcL == 0 && bcR == 0
        plot(x,[0 phi(:,1)' 0],'k')
        axis([min(x) max(x) 0 1.25*max(max(phi))]);
        xlabel('{x(cm)}')
        ylabel('\phi(x)')
    else
        plot(x,phi(:,1),'k')
        axis([min(x) max(x) 0 1.25*max(max(phi))]);
        xlabel('{x(cm)}')
        ylabel('\phi(x)')
    end
end

max(phi(:,1)./phi(:,2))
AL,AD,AU
end

function [k,s,phi,iter] = keffit(c,b,a,dxx,nusig,sct,s,epsK,epsS,itmx )

k       = 1.0;    % guess the initial keff
errK    = 1; 
errS    = 1; 
iter    = 0; 
sindx   = find(s);
n       = length(dxx);
numg    = length(nusig(1,:));
phi     = zeros(n,numg); 

while   ( ( errK>epsK ) || (errS>epsS) ) && iter < itmx
    % call tridiagonal solver --------------------------
    s = s/sum(s);
    phi(:,1)    = tridiag(b(:,1),c(:,1),a(:,1),s/k);
    if numg > 1                     % successive groups
        for g = 2:numg
            scsrc = zeros(n,1);
            for gg = 1:g-1        
            	 scsrc(:,1) = scsrc(:,1) + sct(:,g,gg).*phi(:,gg);
            end 
            % solve for phi_g
            phi(:,g)  = tridiag(c(:,g),b(:,g),a(:,g),scsrc(:,1));
        end
    end
    sold = s; 
    kold = k;   
    s(:) = 0;       % reset fission source
    for g = 1:numg  % compute fission source
      s(:) = s(:) + nusig(:,g).*phi(:,g);
    end  
    % compute keff
    k = sum(s.*dxx)*kold/sum(sold(sindx).*dxx(sindx)); 
    errS = max( abs((s(sindx)-sold(sindx))./s(sindx)) );
    errK = abs( (k-kold)/k );
    iter = iter + 1;  % number of iterations
end
k
disp(['*** final keff estimate: ', num2str(k),' ***'])
disp(['residual errors: errS = ',num2str(errS),' and errK = ',num2str(errK)])
disp(['in ',num2str(iter),' iterations'])

end

